% read in total magnetic field OHCHAMP. 10.02
% convert the julian date & time to local 
% hp filter 11.02
% correlation with mt ts 11.02
% bug & fix history --
%
% latest date 10.12.2002

fid = fopen('\\del00\e\manoj\data\manoj','rt');
% manoj is ascii text, four column ->jd.sec lat log nT
% and very very huge !!!

if fid == -1,
   fprintf('Error in opening a file\n');
   return;
end;

p = 1;
dummy = 1;
seg_size=1024;
   
p=105609204-4000085; % to fix file pointer to some where near 15.7.2001 -> changed 14.6.3
% -from jd 561 and file pointer calculations

fseek(fid,p,-1);

[CA, count] = fscanf(fid,'%[^\n]s',1);
[CA, count] = fscanf(fid,'%f',[4,200000]);
fclose(fid);

% the time axis is jd.as
% jd is the julian date wrt 01.01.2000
% as is normalized seconds ie x/(26*3600)
% ie to get correct seconds multiply as with 26*3600 -> inserted on 09.12.2002
% -could be realized like -> (CA(1,:)-floor(CA(1,:)))/(24*3600)

load \\dell01\e\manoj\projects\champ\kp.dat;

C = length(kp);

for i = 1:C,
   jd(i) = dmy2jd(kp(i,3),kp(i,2),kp(i,1));
end;

jd = jd-dmy2jd(01,01,2000);
C = length(CA);

br(1) = 1;
p=2;

for i = 1:C-1,
if CA(1,i+1)-CA(1,i) > 0.0001,
br(p) = i;
p=p+1;
end;
br(p) = C;
end;


seg_i=0;

[b,a]=cheby2(5,5,0.09,'high');

for i = 1:length(br)-1,
g=CA(4,br(i):br(i+1));
tim=CA(1,br(i):br(i+1));
ll = length(g);

if tim(10)-tim(11)>0, %reverse time axis
g=g(ll:-1:1);
tim=tim(ll:-1:1);
end;

if ll > seg_size,
seg_i=seg_i+1;
crop1 = floor((ll-seg_size)/2); %crop->reason: high lat grad & filter undulations
crop2 = ceil((ll-seg_size)/2);
y=filter(b,a,g);
seg_raw(seg_i,:) = g(crop1+1:ll-crop2);
seg_flt(seg_i,:) = y(crop1+1:ll-crop2);
seg_tim(seg_i,:) = tim(crop1+1:ll-crop2);
end;
end;
%return;

kk=[9 10 11 13 14 16 17 18 19 21 22 23];

for i = 1:length(br),
   T = floor(seg_tim(i,1));
   T1 = kp(find(jd==T),4);
   subplot(4,4,find(kk==T1));
   %d=std(seg_flt(i,:));
   plot(seg_flt(i,:));
   
   %plot(T1,d,'r.');
   axis([-inf inf -inf inf]);
	ylabel('nT');
   xlabel(['Kp -' sprintf('%d',T1) 'JD -' sprintf('%d',T)]);
   hold on;
    
end;